######################################################
##Replication file for "Hafner-Burton, Hyde and Jablonski. 
##Surviving Elections: Election Violence, Incumbent Victory, and Post-Election Repercussions" 
##British Journal of Political Science. Accepted November 2015
##
##Plots Figure 4
##Run on R 2.14.2
######################################################

rm(list = ls())
library(foreign)
library(Zelig)
library(gee)
library(ggplot2)
library(foreign)
library(nlme)
library(sandwich)
library(car)
require(grid)
require(lmtest)
library(boot)
library(AICcmodavg)
library(meboot)
set.seed(1234)

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

tenure.data<-read.dta("Stata12BJPSReplication12112015.dta")

tenure.data<-subset(tenure.data, !is.na(tenure.data$vt))
tenure.data<-subset(tenure.data, !is.na(tenure.data$uncertain))
tenure.data<-subset(tenure.data, !is.na(tenure.data$otherpreviolence))
tenure.data<-subset(tenure.data, !is.na(tenure.data$ma_physint))
tenure.data<-subset(tenure.data, !is.na(tenure.data$age))
tenure.data<-subset(tenure.data, !is.na(tenure.data$sumten))
tenure.data<-subset(tenure.data, !is.na(tenure.data$civwar))
tenure.data<-subset(tenure.data, !is.na(tenure.data$gdp))
tenure.data<-subset(tenure.data, !is.na(tenure.data$pop))
tenure.data<-subset(tenure.data, !is.na(tenure.data$ma_polity2))
tenure.data<-subset(tenure.data, !is.na(tenure.data$any_demonstration))
tenure.data<-subset(tenure.data, !is.na(tenure.data$prefraud))
tenure.data<-subset(tenure.data, !is.na(tenure.data$compulsory))
tenure.data<-subset(tenure.data, !is.na(tenure.data$multipleround))


attach(tenure.data)
tenure.data.incumbents = subset(tenure.data, subset=(incumbentrun==1 & competitive==1 & finalround==1))
detach(tenure.data)

turnout.range = seq(min(tenure.data.incumbents$vt), max(tenure.data.incumbents$vt), by=.1)


z.out<-zelig(formula=incumbentwins~vt + turnoutviolence + otherpreviolence + ma_physint + age + sumten + civwar + gdp + pop + uncertain + ma_polity2 + any_demonstration + prefraud + compulsory + multipleround, data=tenure.data.incumbents, model="logit.gee",  id="ccode")
summary(z.out)


nsims=1000
yhat0=matrix(nrow=nsims, ncol=length(turnout.range))
for (i in 1:length(turnout.range) ) {
	x.out.range <- setx(z.out, turnoutviolence=turnout.range[i]*0, otherpreviolence=0, vt=turnout.range[i])
	s.out.range <- sim(z.out, x=x.out.range, num=nsims)
	yhat0[,i] = s.out.range$qi$ev1
}

yhat1=matrix(nrow=nsims, ncol=length(turnout.range))
for (i in 1:length(turnout.range) ) {
	x.out.range <- setx(z.out, turnoutviolence=turnout.range[i]*1, otherpreviolence=1, vt=turnout.range[i])
	s.out.range <- sim(z.out, x=x.out.range, num=nsims)
	yhat1[,i] = s.out.range$qi$ev1
}

predict.data0=data.frame(cbind(upper=apply(yhat0, 2, quantile, probs=.95), lower=apply(yhat0, 2, quantile, probs=.05), yhat=apply(yhat0, 2, mean), treat=seq(min(tenure.data.incumbents$vt), max(tenure.data.incumbents$vt), by=.1), sd=apply(yhat0, 2, sd)))
predict.data1=data.frame(cbind(upper=apply(yhat1, 2, quantile, probs=.95), lower=apply(yhat1, 2, quantile, probs=.05), yhat=apply(yhat1, 2, mean), treat=seq(min(tenure.data.incumbents$vt), max(tenure.data.incumbents$vt), by=.1), sd=apply(yhat1, 2, sd)))

predictplot0=ggplot(predict.data0, aes(x=treat, y=yhat)) + 
geom_linerange(ymax=predict.data0$upper, ymin=predict.data0$lower, colour="grey") +
geom_line(colour="black", width=2, size=2) + 
theme_bw(base_family="serif", base_size=20) +
xlab("Voter Turnout") +
ylab("Pr(Incumbent Wins)") +
scale_y_continuous(limits=c(0, 1)) +
scale_x_continuous(breaks=c(0, 100), labels=c("MIN\n 0", "MAX\n 100")) +
theme(
	panel.grid.major=element_blank(),
	panel.grid.minor=element_blank(), 
	panel.border = element_blank(), 
	panel.background = element_blank(),
	axis.line = element_line(colour = "black", size=1),
	plot.margin =  unit(c(1,1,1,1), "cm"),
	axis.title.x = element_text(vjust = 2),
	axis.title.y = element_text(vjust = .25),
	axis.text.x = element_text(size=19,vjust = 0),
	axis.ticks.x = element_line(colour = "black", size=1),
	axis.ticks.length = unit(c(.2), "cm")
	) 

predictplot0=predictplot0 + ggtitle("Pre-Election Violence=0")


predictplot1=ggplot(predict.data1, aes(x=treat, y=yhat)) + 
geom_linerange(ymax=predict.data1$upper, ymin=predict.data1$lower, colour="grey") +
geom_line(colour="black", width=2, size=2) + 
theme_bw(base_family="serif", base_size=20) +
xlab("Voter Turnout") +
ylab("Pr(Incumbent Wins)") +
scale_y_continuous(limits=c(0, 1)) +
scale_x_continuous(breaks=c(0, 100), labels=c("MIN\n 0", "MAX\n 100")) +
theme(
	panel.grid.major=element_blank(),
	panel.grid.minor=element_blank(), 
	panel.border = element_blank(), 
	panel.background = element_blank(),
	axis.line = element_line(colour = "black", size=1),
	plot.margin =  unit(c(1,1,1,1), "cm"),
	axis.title.x = element_text(vjust = 2),
	axis.title.y = element_text(vjust = .25),
	axis.text.x = element_text(size=19,vjust = 0),
	axis.ticks.x = element_line(colour = "black", size=1),
	axis.ticks.length = unit(c(.2), "cm")
	) 

predictplot1=predictplot1 + ggtitle("Pre-Election Violence=1")

multiplot(predictplot1, predictplot0, cols=2)